
#include "errlog.h"
#include "grid3d.h"
// #include "alegrid.h"
#include "f3_parser.h"
#include "glovar.h"
#include "fe.h"
#include "fe_matr.h"
#include <string.h>
#include <stdlib.h>

/************ The Flags *******************************/
struct f3_flags_t flags;
const unsigned per_default=per_broyden;
const char * const periodic_str[6]={"off", "marching", "adjoint", "armijo", "broyden", "lbfgs"};

/************ The Physical Parameters *****************/
double Re;
double Fr;

void get_coefs(unsigned advection, double *mass_coef, double *stif_coef, double *advm_coef)
{
    if(mass_coef) *mass_coef=1.0;
    if(stif_coef) *stif_coef=TStep;
    if(advm_coef)
        switch(advection)
        {
            case adv_none: *advm_coef=0.0;
                break; /* no advection (Stokes) */
            case adv_explicit: *advm_coef=0.0;
                break; /* explicit advection */
            case adv_semi: *advm_coef=TStep;
                break; /* semi-implicit advection */
            case adv_char: *advm_coef=0.0;
                break; /* characteristics advection */
            case adv_rk45: *advm_coef=0.0;
                break; /* runge-kutta */
        }
    return;
}

/************ The Time Parameters **********************/
int total_ti, curr_ti;
double TStep, StopTime;

void SetTimeParams(double start, double stop, double step)
{
    PrintLog("glovar: SetTimeParams: ");
    if((start<0.0)||(stop<=start))
        error_msg("SetTimeParams: invalid start and/or stop times.\n");
    if(step<=0||step>stop-start)
    {
        step=stop-start;
        warning_msg("SetTimeParams: step must be positive and <= stop-stat. fixed.\n");
    }
    StopTime=stop;
    TStep=step;
    total_ti=floor(StopTime/TStep);
    if(StopTime-total_ti*TStep>SMALLF) ++total_ti;
    curr_ti=floor(start/TStep);
    PrintLog("TStep=%-10f  StopTime=%-10f\n", (double) TStep, (double) StopTime);
}

/************ The Grid *********************************/
// alegrid_t *alegrid=NULL;
grid3d_t *grid=NULL;
mesq_handle mesquite;
node_t *adv_nodes=NULL;

/************ The solution vectors *********************/
vector u, u0, u1, adv_u;

scalar p, p0, p1, phi;

/************ The adjoint solution *********************/
vector v, v0, v1;

scalar q, q0, q1, qhi;

/************ some matrices ****************************/
fe_matr_p massP=NULL,
        massU=NULL;

void alloc_glovar(unsigned elemU, unsigned elemP)
{
    int vec_dim, scal_dim;

    scal_dim=fe_ndof(grid, elemP);
    if(elemU==fe_bub_vec)
        vec_dim=fe_ndof(grid, fe_bub);
    else
        vec_dim=fe_ndof(grid, elemU);

    u=alloc_vec(vec_dim);
    u0=alloc_vec(vec_dim);
    u1=alloc_vec(vec_dim);
    adv_u=alloc_vec(vec_dim);

    p=alloc_scal(scal_dim);
    p0=alloc_scal(scal_dim);
    p1=alloc_scal(scal_dim);
    phi=alloc_scal(scal_dim);

    v=alloc_vec(vec_dim);
    v0=alloc_vec(vec_dim);
    v1=alloc_vec(vec_dim);

    q=alloc_scal(scal_dim);
    q0=alloc_scal(scal_dim);
    q1=alloc_scal(scal_dim);
    qhi=alloc_scal(scal_dim);

    massP=fe_matr_alloc(elemP, local_mass, eo_add_all, grid, TRUE,
#if defined(USE_CSC) 
            "ss"
#else
            NULL
#endif
            );
    // fe_matr_structure(grid, 1, &massP);
    fe_matr_assemble(grid, 1, &massP);

    massU=fe_matr_alloc(elemU, local_mass, eo_add_all, grid, TRUE,
#if defined(USE_CSC) 
            "ss"
#else
            NULL
#endif
            );
    // fe_matr_structure(grid, 1, &massU);
    fe_matr_assemble(grid, 1, &massU);

}

void clean_glovecs(void)
{
    clean_vec(u);
    clean_vec(u0);
    clean_vec(u1);
    clean_vec(adv_u);

    clean_scal(p);
    clean_scal(p0);
    clean_scal(p1);
    clean_scal(phi);

    clean_vec(v);
    clean_vec(v0);
    clean_vec(v1);

    clean_scal(q);
    clean_scal(q0);
    clean_scal(q1);
    clean_scal(qhi);
}

void free_glovar(void)
{
    free_vec(u);
    free_vec(u0);
    free_vec(u1);
    free_vec(adv_u);
    free_scal(phi);
    free_scal(p);
    free_scal(p0);
    free_scal(p1);

    free_vec(v);
    free_vec(v0);
    free_vec(v1);
    free_scal(q);
    free_scal(q0);
    free_scal(q1);
    free_scal(qhi);

    fe_matr_free(massP);
    fe_matr_free(massU);
    if(adv_nodes)
        free(adv_nodes);
}


/************ The Functions ****************************/

void (*exct_vel)(double, const double *, double *)=NULL;
void (*exct_prs)(double, const double *, double *)=NULL;
void (*init_vel)(double, const double *, double *)=NULL;
void (*init_prs)(double, const double *, double *)=NULL;
void (*srce_trm)(double, const double *, double *)=NULL;

void (*adj_exct_vel)(double, const double *, double *)=NULL;
void (*adj_exct_prs)(double, const double *, double *)=NULL;
void (*adj_finl_vel)(double, const double *, double *)=NULL;
void (*adj_finl_prs)(double, const double *, double *)=NULL;
void (*adj_srce_trm)(double, const double *, double *)=NULL;

static char *func_strings[10];

enum
{
    e_vel_ind=0,
    e_prs_ind=1,
    i_vel_ind=2,
    i_prs_ind=3,
    s_trm_ind=4,
    adj_e_vel_ind=5,
    adj_e_prs_ind=6,
    adj_f_vel_ind=7,
    adj_f_prs_ind=8,
    adj_s_trm_ind=9
};

static void _glovar_parsing_error(const char *func_name)
{
    error_msg("Failed to parse %s string.\n", func_name);
}

static void _e_vel(double T, const double *Pt, double *Val)
{
    if(func_strings[e_vel_ind]==NULL)
        Val[0]=Val[1]=Val[2]=0.0;
    else if(!parse_vect(T, Pt, Val, func_strings[e_vel_ind]))
        _glovar_parsing_error("exct_vel");
}

static void _e_prs(double T, const double *Pt, double *Val)
{
    if(func_strings[e_prs_ind]==NULL)
        Val[0]=0.0;
    else if(!parse_scal(T, Pt, Val, func_strings[e_prs_ind]))
        _glovar_parsing_error("exct_prs");
}

static void _i_vel(double T, const double *Pt, double *Val)
{
    if(func_strings[i_vel_ind]==NULL)
        Val[0]=Val[1]=Val[2]=0.0;
    else if(!parse_vect(T, Pt, Val, func_strings[i_vel_ind]))
        _glovar_parsing_error("init_vel");
}

static void _i_prs(double T, const double *Pt, double *Val)
{
    if(func_strings[i_prs_ind]==NULL)
        Val[0]=0.0;
    else if(!parse_scal(T, Pt, Val, func_strings[i_prs_ind]))
        _glovar_parsing_error("init_prs");
}

static void _s_trm(double T, const double *Pt, double *Val)
{
    if(func_strings[s_trm_ind]==NULL)
        Val[0]=Val[1]=Val[2]=0.0;
    else if(!parse_vect(T, Pt, Val, func_strings[s_trm_ind]))
        _glovar_parsing_error("srce_trm");
}

static void _adj_e_vel(double T, const double *Pt, double *Val)
{
    if(func_strings[adj_e_vel_ind]==NULL)
        Val[0]=Val[1]=Val[2]=0.0;
    else if(!parse_vect(T, Pt, Val, func_strings[adj_e_vel_ind]))
        _glovar_parsing_error("adj_exct_vel");
}

static void _adj_e_prs(double T, const double *Pt, double *Val)
{
    if(func_strings[adj_e_prs_ind]==NULL)
        Val[0]=0.0;
    else if(!parse_scal(T, Pt, Val, func_strings[adj_e_prs_ind]))
        _glovar_parsing_error("adj_exct_prs");
}

static void _adj_f_vel(double T, const double *Pt, double *Val)
{
    if(func_strings[adj_f_vel_ind]==NULL)
        Val[0]=Val[1]=Val[2]=0.0;
    else if(!parse_vect(T, Pt, Val, func_strings[adj_f_vel_ind]))
        _glovar_parsing_error("adj_finl_vel");
}

static void _adj_f_prs(double T, const double *Pt, double *Val)
{
    if(func_strings[adj_f_prs_ind]==NULL)
        Val[0]=0.0;
    else if(!parse_scal(T, Pt, Val, func_strings[adj_f_prs_ind]))
        _glovar_parsing_error("adj_finl_prs");
}

static void _adj_s_trm(double T, const double *Pt, double *Val)
{
    if(func_strings[adj_s_trm_ind]==NULL)
        Val[0]=Val[1]=Val[2]=0.0;
    else if(!parse_vect(T, Pt, Val, func_strings[adj_s_trm_ind]))
        _glovar_parsing_error("adj_srce_trm");
}

static int _set_func_str(int indx, const char *f_str)
{
    char *s;
    int len;
    if(f_str)
    {
        if(f_str[0]=='\0')
        {
            if(func_strings[indx])
            {
                free(func_strings[indx]);
                func_strings[indx]=NULL;
            }
            return FALSE;
        }
        len=strlen(f_str)+1;
        s=(char*) realloc(func_strings[indx], len*sizeof(char));
        if(s==NULL)
            error_msg("Memory allocation error in SetFuncStrings.\n");
        func_strings[indx]=strncpy(s, f_str, len);
        return TRUE;
    }
    else
        return func_strings[indx]!=NULL;
}

void SetFuncStrings(const char *init_vel_str,
                    const char *init_prs_str, const char *srce_trm_str)
{
    init_vel=_set_func_str(i_vel_ind, init_vel_str)?_i_vel:NULL;
    init_prs=_set_func_str(i_prs_ind, init_prs_str)?_i_prs:NULL;
    srce_trm=_set_func_str(s_trm_ind, srce_trm_str)?_s_trm:NULL;
}

void SetAdjFuncStrings(const char *adj_finl_vel_str,
                       const char *adj_finl_prs_str, const char *adj_srce_trm_str)
{
    adj_finl_vel=_set_func_str(adj_f_vel_ind, adj_finl_vel_str)?_adj_f_vel:NULL;
    adj_finl_prs=_set_func_str(adj_f_prs_ind, adj_finl_prs_str)?_adj_f_prs:NULL;
    adj_srce_trm=_set_func_str(adj_s_trm_ind, adj_srce_trm_str)?_adj_s_trm:NULL;
}

void SetExctStrings(const char *exct_vel_str, const char *exct_prs_str)
{
    exct_vel=_set_func_str(e_vel_ind, exct_vel_str)?_e_vel:NULL;
    exct_prs=_set_func_str(e_prs_ind, exct_prs_str)?_e_prs:NULL;
}

void SetAdjExctStrings(const char *adj_exct_vel_str, const char *adj_exct_prs_str)
{
    adj_exct_vel=_set_func_str(adj_e_vel_ind, adj_exct_vel_str)?_adj_e_vel:NULL;
    adj_exct_prs=_set_func_str(adj_e_prs_ind, adj_exct_prs_str)?_adj_e_prs:NULL;
}

/****************************************************************************
 * Setting boundary conditions
 */

typedef struct bvd_s
{
    struct bvd_s *next;
    unsigned on_bdry;
    char *f_string;
    field3d *f_routine;
} bvd_t;

void (*bdry_vel)(unsigned, double, const double *, double *)=NULL;
void (*adj_bdry_vel)(unsigned, double, const double *, double *)=NULL;

static void _b_vel_generic(bvd_t *root, unsigned on_bdry,
                           double T, const double *Pt, double *V)
{
    bvd_t *bvd;

    for(bvd=root; bvd!=NULL; bvd=bvd->next)
    {
        if(bvd->on_bdry&on_bdry)
        {
            if(bvd->f_routine)
                bvd->f_routine(T, Pt, V);
            else if(bvd->f_string)
            {
                if(!parse_vect(T, Pt, V, bvd->f_string))
                    _glovar_parsing_error("bdry_vel");
            }
            else
                V[0]=V[1]=V[2]=0.0;
            return;
        }
    }
    V[0]=V[1]=V[2]=0.0;
    return;
}

static bvd_t *_get_bvd_genedic(bvd_t **root, unsigned on_bdry,
                               void (**vel)(unsigned, double, const double *, double *),
                               void (*_vel)(unsigned, double, const double *, double *)
                               )
{
    if(*root==NULL)
    {
        *root=calloc(1, sizeof(bvd_t));
        if(*root==NULL)
            error_msg("Memory allocation error in %s\n", __func__);
        (*root)->on_bdry=on_bdry;
        *vel=_vel;
        return *root;
    }
    else
    {
        bvd_t *bvd= *root;
        while(bvd->next)
        {
            if(bvd->on_bdry==on_bdry)
                return bvd;
            bvd=bvd->next;
        }
        if(bvd->on_bdry==on_bdry)
            return bvd;
        bvd->next=calloc(1, sizeof(bvd_t));
        if(bvd->next==NULL)
            error_msg("Memory allocation error in %s\n", __func__);
        bvd=bvd->next;
        bvd->on_bdry=on_bdry;
        return bvd;
    }
}

static bvd_t *bvd_root=NULL;

static void _b_vel(unsigned on_bdry, double T, const double *Pt, double *V)
{
    _b_vel_generic(bvd_root, on_bdry, T, Pt, V);
}

static bvd_t *_get_bvd(unsigned on_bdry)
{
    return _get_bvd_genedic(&bvd_root, on_bdry, &bdry_vel, _b_vel);
}

static bvd_t *adj_bvd_root=NULL;

static void _adj_b_vel(unsigned on_bdry, double T, const double *Pt, double *V)
{
    _b_vel_generic(adj_bvd_root, on_bdry, T, Pt, V);
}

static bvd_t *_get_adj_bvd(unsigned on_bdry)
{
    return _get_bvd_genedic(&adj_bvd_root, on_bdry, &adj_bdry_vel, _adj_b_vel);
}

void SetBdryFuncString(unsigned on_bdry, const char *bv_str)
{
    int len;
    char *s;
    bvd_t *bvd;

    bvd=_get_bvd(on_bdry);
    bvd->f_routine=NULL;
    if(bvd->f_string)
        free(bvd->f_string);
    if(bv_str==NULL||bv_str[0]=='\0')
    {
        bvd->f_string=NULL;
        return;
    }
    len=strlen(bv_str)+1;
    s=(char *) calloc(len, sizeof(char));
    if(s==NULL)
        error_msg("Memory allocation error in %s.\n", __func__);
    bvd->f_string=strncpy(s, bv_str, len);
    return;
}

void SetAdjBdryFuncString(unsigned on_bdry, const char *bv_str)
{
    int len;
    char *s;
    bvd_t *bvd;

    bvd=_get_adj_bvd(on_bdry);
    bvd->f_routine=NULL;
    if(bvd->f_string)
        free(bvd->f_string);
    if(bv_str==NULL||bv_str[0]=='\0')
    {
        bvd->f_string=NULL;
        return;
    }
    len=strlen(bv_str)+1;
    s=(char *) calloc(len, sizeof(char));
    if(s==NULL)
        error_msg("Memory allocation error in %s.\n", __func__);
    bvd->f_string=strncpy(s, bv_str, len);
    return;
}

void SetBdryFuncRoutine(unsigned on_bdry, field3d *b_vel)
{
    bvd_t *bvd;

    bvd=_get_bvd(on_bdry);
    bvd->f_routine=b_vel;
    if(bvd->f_string)
    {
        free(bvd->f_string);
        bvd->f_string=NULL;
    }
    return;
}

void SetAdjBdryFuncRoutine(unsigned on_bdry, field3d *adj_b_vel)
{
    bvd_t *bvd;

    bvd=_get_adj_bvd(on_bdry);
    bvd->f_routine=adj_b_vel;
    if(bvd->f_string)
    {
        free(bvd->f_string);
        bvd->f_string=NULL;
    }
    return;
}


/*****************************************************************************/

void (*shgrd_vel)(unsigned, double, const double *, double *)=NULL;
static bvd_t *shgrd_root=NULL;

static void _shgrd_vel(unsigned on_bdry, double T, const double *Pt, double *V)
{
    _b_vel_generic(shgrd_root, on_bdry, T, Pt, V);
}

static bvd_t *_get_shgrd_bvd(unsigned on_bdry)
{
    return _get_bvd_genedic(&shgrd_root, on_bdry, &shgrd_vel, _shgrd_vel);
}

void SetShGrdBdryFuncString(unsigned on_bdry, const char *sg_str)
{
    int len;
    char *s;
    bvd_t *bvd;

    bvd=_get_shgrd_bvd(on_bdry);
    bvd->f_routine=NULL;
    if(bvd->f_string)
        free(bvd->f_string);
    if(sg_str==NULL||sg_str[0]=='\0')
    {
        bvd->f_string=NULL;
        return;
    }
    len=strlen(sg_str)+1;
    s=(char *) calloc(len, sizeof(char));
    if(s==NULL)
        error_msg("Memory allocation error in %s.\n", __func__);
    bvd->f_string=strncpy(s, sg_str, len);
    return;
}

void SetShGrdBdryFuncRoutine(unsigned on_bdry, field3d *sg_vel)
{
    bvd_t *bvd;

    bvd=_get_shgrd_bvd(on_bdry);
    bvd->f_routine=sg_vel;
    if(bvd->f_string)
    {
        free(bvd->f_string);
        bvd->f_string=NULL;
    }
    return;
}
